/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package ProteoXL;

import java.io.*;
import java.text.DecimalFormat;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.Iterator;
import java.util.List;
import uk.ac.ebi.jmzml.model.mzml.Spectrum;
import uk.ac.ebi.jmzml.model.mzml.*;
import uk.ac.ebi.jmzml.xml.io.MzMLObjectIterator;
import uk.ac.ebi.jmzml.xml.io.MzMLUnmarshaller;

/**
 * This class will read two inputs
 * - detected feature masses e.g. from Progenesis, TransOmics or else
 * - Output of DeconvTools to generate a peak list based on those MSe ions sharing the same RT profile
 *      - We need to know what is MS1 and what is MSn/MSe -- may be better to run DeconvTools on function=2 spectra only
 *      - Perl script cannot be located that does this work - recode using mzML directly here to access metadata?
 * 
 * - This is going to produce an MGF file for searching - all fragments across scans that are within the bounds of an MS1 feature, are assigned as putatively belonging to it
 * 
 * To Do:
 * 
 * --> Removal of noise -- currently all peaks included
 * --> May want to improve matching of fragments to precursors based on chromatogram shape - currently just using APEX
 *      --> This method may bias in favour of noise over genuine peaks, since noise peaks will only appear once
 *      --> May want to have seen fragment more than once over rt window to allow it as a fragment
 *      --> May also need to optimise the running of DeconvTools
 * 
 * - 
 * @author jonesar
 */
public class ProcessDeconv {

    private static final double PROTONMASS = 1.00727647;
    String mzMLFile = "C:/Work/Crosslinking/Waters/20110505_HW_BSATryptic_SS_mapping_03.mzML";
    //String mzMLFile = "C:/Work/Crosslinking/Waters/BSA_subset/First1000seconds/20110505_HW_BSATryptic_SS_mapping_03.mzML";
    String mgfFile = "C:/Work/Crosslinking/Waters/DeconvToolsOutput/20110505_HW_BSATryptic_SS_mapping_03_decharged.mgf";
    //String mgfFile = "C:/Work/Crosslinking/Waters/BSA_subset/First1000seconds/20110505_HW_BSATryptic_SS_mapping_03_decharged.mgf";
    String deconvToolsScanTable = "C:/Work/Crosslinking/Waters/DeconvToolsOutput/20110505_HW_BSATryptic_SS_mapping_03_scans.csv";
    
    DecimalFormat df = new DecimalFormat("#.#");
    
    private static int RT_WIN_FOR_MATCHING_FRAG_APEX = 3;   //i.e. APEX of fragments has to be within 3 seconds of feature APEX
    
        /*
     * 
     * scan_num	scan_time	type	bpi	bpi_mz	tic	num_peaks	num_deisotoped	info
        0	126.516	1	1145	785.84362	17747	43	2	0
        1	127.565	1	995	153.13822	30530	141	0	1

     */
    
    String readMS1FeaturesFile = "C:/Work/Crosslinking/Waters/BSA_TransOmics/features_transomics.csv";
    /*
     * 								Normalized abundance	Raw abundance	Intensity	Sample retention time (min)	
								Condition 1	Condition 1	Condition 1	Condition 1	
#	m/z	Retention time (min)	Retention time window (min)	Rt seconds	Mass	Charge	Maximum CV	20110505_HW_BSATryptic_SS_mapping_03	20110505_HW_BSATryptic_SS_mapping_03	20110505_HW_BSATryptic_SS_mapping_03	20110505_HW_BSATryptic_SS_mapping_03	Notes
445	674.2707577	3.232366667	0.122333333	193.942	1346.526962	2		296.4256525	296.4256525	1261	3.232366667	
313	390.2151857	3.907783333	0.122333333	234.467	389.2079093	1		506.8514848	506.8514848	2922	3.907783333	
178	537.2843674	3.907783333	0.174766667	234.467	536.277091	1		1232.391903	1232.391903	5311	3.907783333	
239	717.3233548	5.566966667	0.343333333	334.018	716.3160783	1		2300.398022	2300.398022	3772	5.566966667	

     * 
     */
    
    
    ArrayList<Object[]> ms1Features = new ArrayList();  //Object[] is mass, charge, rt and rt window
    String deconvToolsIsoTable = "C:/Work/Crosslinking/Waters/DeconvToolsOutput/20110505_HW_BSATryptic_SS_mapping_03_isos.csv";
    
    /*
     * scan_num	charge	abundance	mz	fit	average_mw	monoisotopic_mw	mostabundant_mw	fwhm	signal_noise	mono_abundance	mono_plus2_abundance	flag	interference_score
        0	2	1145	785.84362	0.0888	1570.72075	1569.67378	1569.67378	0.0326	42.42	1145	365		0
        0	2	155	777.34023	0.1973	1552.68858	1551.65363	1551.65363	0.0242	7.34	155	66		0
        2	1	1035	153.1401	0.0125	152.22856	152.13338	152.13338	0.0119	57.22	1035	0		0
* 
     * 
     */
    
    HashMap<String,Spectrum> rt1dpToSpectra = new HashMap();
    HashMap<Integer,String> scanNumToRT1DP = new HashMap();
    HashMap<String,Integer> rt1dpToScanNum = new HashMap();
    HashMap<Integer,ArrayList<Object[]>> scanNumToAllDechargedData = new HashMap(); //Scan number references an arraylist of data points - each data point is mass and intensity pairs
    
    HashMap<Integer,Integer> rtIntegerToMSNScanNum = new HashMap();    //TODO - this works off the assumption that only one MSn scan will have a given RT (integer value)
    
    
    
       /**
     * @param args the command line arguments
     */
    public static void main(String[] args) {
        
        ProcessDeconv deconv = new ProcessDeconv();
        
    }
    
    
    public ProcessDeconv() {
        
        readFeaturesFile(readMS1FeaturesFile);
        readDeconvScanTable(deconvToolsScanTable);
        readDeconvIsoTable(deconvToolsIsoTable);
        readSpectra(mzMLFile);
        GenerateDechargedSpectrum(mgfFile);
    }
    
    

    public void readSpectra(String xmlFile){
        final String xmlfile = xmlFile;
        File file = new File(xmlfile);
        if (file.exists())
        {                
            MzMLUnmarshaller unmarshaller = new MzMLUnmarshaller(file);
            MzMLObjectIterator<Spectrum> spectrumIterator = unmarshaller.unmarshalCollectionFromXpath("/run/spectrumList/spectrum", Spectrum.class);
            while (spectrumIterator.hasNext()){
                Spectrum spectrum = spectrumIterator.next();
                String mslevel = "";                        
                List<CVParam> specParam = spectrum.getCvParam();      
                for (Iterator lCVParamIterator = specParam.iterator(); lCVParamIterator.hasNext();){
                    CVParam lCVParam = (CVParam) lCVParamIterator.next();
                    if (lCVParam.getAccession().equals("MS:1000511")){
                        mslevel = lCVParam.getValue().trim();
                    }
                }
                double rt = 0.0;
                String unitRT = "";
                List<CVParam> scanParam = spectrum.getScanList().getScan().get(0).getCvParam();
                for (Iterator lCVParamIterator = scanParam.iterator(); lCVParamIterator.hasNext();){
                    CVParam lCVParam = (CVParam) lCVParamIterator.next();
                    if (lCVParam.getAccession().equals("MS:1000016")){
                        unitRT = lCVParam.getUnitAccession().trim();
                        if (unitRT.equals("UO:0000031")) //... Validating rt unit (mins or secs) ...//
                        {    
                            rt = Float.parseFloat(lCVParam.getValue().trim()) * 60;
                        }
                        else{
                            rt = Float.parseFloat(lCVParam.getValue().trim());
                        }
                    }  
                }
                //model.insertRow(model.getRowCount(), new Object[]{spectrum.getIndex(),spectrum.getId(),mslevel,Float.parseFloat(String.format("%.2f", rt))});
                
                String specID = spectrum.getId();
                String formatRT = String.format("%.2f", rt);                    
                rt1dpToSpectra.put(formatRT, spectrum);
                Integer scanNum = rt1dpToScanNum.get(formatRT);
                
                if(specID.indexOf("function=2")!=-1){   //Get function=2 spectra only

                    //System.out.println("MSn\tspecID:" + spectrum.getId() + " level: " + mslevel + " rt: " + formatRT + "scan num: " + scanNum);
                    //Now let's get the peaks for this spectrum
                    
                    if(scanNum != null){
                        
                        ArrayList<Object[]> peaks = scanNumToAllDechargedData.get(scanNum);
                        
                        /*
                        if(peaks!= null){
                            for(Object[] peak : peaks){
                                double mass = (double)peak[0];
                                int intensity = (int)peak[1];
                                System.out.println("\t" + mass + " " + intensity);
                            }
                        }
                        */
                        rtIntegerToMSNScanNum.put(new Integer((int)rt), scanNum);   //Cast rt to integer for later retrieval
                    }
                }
                else{
                    //System.out.println("MS1\tspecID:" + spectrum.getId() + " level: " + mslevel + " rt: " + formatRT + "scan num: " + scanNum);
                    
                    
                }
                        //function=2
            }

        }
        else
        {
            System.out.println("File not found");
        }
                                        
    } 
    
    /*
     * Read the scan table produced by deconv tools
     *      * scan_num	scan_time	type	bpi	bpi_mz	tic	num_peaks	num_deisotoped	info
        0	126.516	1	1145	785.84362	17747	43	2	0
        1	127.565	1	995	153.13822	30530	141	0	1
          * 
     */
    public void readDeconvScanTable(String infile){

        BufferedReader br = null;
 
        try {
            String line = null;
            br = new BufferedReader(new FileReader(infile));

            int i = 0;
            while ((line = br.readLine()) != null) {
                    //System.out.println(sCurrentLine);
                
                if(i!=0){//ignore header line
                    String[] temp = line.split(",");
                    Integer scanNum = Integer.parseInt(temp[0]);
                    String rtString = String.format("%.2f", Float.parseFloat(temp[1])); //Format to 2dp
                    //System.out.println("ScanNum:" + scanNum + " rt: " + rtString);
                    scanNumToRT1DP.put(scanNum, rtString);
                    rt1dpToScanNum.put(rtString,scanNum);
                }
                i++;
                    
            }
        } 
        catch (IOException e) {
                e.printStackTrace();
        } 
        
    }
    
        
    
    /*
     * Read table in this kind of format and associate all rows of data with right scan number
     * scan_num	charge	abundance	mz	fit	average_mw	monoisotopic_mw	mostabundant_mw	fwhm	signal_noise	mono_abundance	mono_plus2_abundance	flag	interference_score
        0	2	1145	785.84362	0.0888	1570.72075	1569.67378	1569.67378	0.0326	42.42	1145	365		0
        0	2	155	777.34023	0.1973	1552.68858	1551.65363	1551.65363	0.0242	7.34	155	66		0
        2	1	1035	153.1401	0.0125	152.22856	152.13338	152.13338	0.0119	57.22	1035	0		0
    * 
     * 
     */
    public void readDeconvIsoTable(String infile){

        BufferedReader br = null; 
        try {
            String line = null;
            br = new BufferedReader(new FileReader(infile));

            int i = 0;
            while ((line = br.readLine()) != null) {
                
                if(i!=0){//ignore header line
                    String[] temp = line.split(",");
                    Integer scanNum = Integer.parseInt(temp[0]);
                    Integer abundance = Integer.parseInt(temp[2]);      //This is designed for Waters data only - i.e. integer ion counts
                    int charge = Integer.parseInt(temp[1]);
                    Double mz = Double.parseDouble(temp[3]);
                    Double mass = (mz * charge) - (charge * PROTONMASS);
                    
                    Object[] peak = {mass,abundance};
                    
                    ArrayList<Object[]> peaksForThisSpectrum = null;
                    if(scanNumToAllDechargedData.containsKey(scanNum)){
                        peaksForThisSpectrum = scanNumToAllDechargedData.get(scanNum);
                    }
                    else{
                        peaksForThisSpectrum = new ArrayList();
                        scanNumToAllDechargedData.put(scanNum, peaksForThisSpectrum);
                    }
                    peaksForThisSpectrum.add(peak);
                    
                }
                i++;
                    
            }
        } 
        catch (IOException e) {
                e.printStackTrace();
        }         
    }
    
    
    
    //#	m/z	Retention time (min)	Retention time window (min)	Rt seconds	Mass
    public void readFeaturesFile(String infile){
        
        BufferedReader br = null; 
        try {
            String line = null;
            br = new BufferedReader(new FileReader(infile));

            HashMap<String,Integer> headersToCols = new HashMap();//map column headers to columns
            int i = 0;
            while ((line = br.readLine()) != null) {
                String[] temp = line.split(",");
                
                if(i==2){
                    int j=0;
                    for(String header : temp){
                        headersToCols.put(temp[j], new Integer(j));
                        j++;
                    }
                }
                if(i>=3){//ignore 3header lines before data
                    
                    Double mass = Double.parseDouble(temp[headersToCols.get("Mass")]);
                    int charge = Integer.parseInt(temp[headersToCols.get("Charge")]);
                    Double rt = Double.parseDouble(temp[headersToCols.get("Retention time (min)")])*60;
                    Double rtWindow = Double.parseDouble(temp[headersToCols.get("Retention time window (min)")])*60;
                    
                    Object[] feature = {mass,charge,rt,rtWindow};
                    ms1Features.add(feature);      
                    //System.out.println("Feature" + mass + " " + charge + " " + rt + " " + rtWindow);
                    
                }
                i++;
                    
            }
        } 
        catch (IOException e) {
                e.printStackTrace();
        } 
        
    }
    
    /*
     * Method will do following:
     * - For each found MS1 feature:
     * - Get RT window (presume window is half in either direction
     * - Extract all MSn peaks with mass < feature mass, and create a spectrum
     */
    private void GenerateDechargedSpectrum(String outFile){
        
        try{
            FileWriter writer = new FileWriter(outFile);
            
            for(Object[] feature : ms1Features){
                //rt1dpToScanNum - try various lookups into the hashmap to get the spectra by integer
                //Object[] feature = {mass,charge,rt,rtWindow};

                double mass = (double)feature[0];
                int charge = (int)feature[1];
                double peptideRT = (double)feature[2];
                double peptideRTWin = (double)feature[3];

                int rtStart = (int)(peptideRT - (peptideRTWin/2));
                int rtEnd = (int)(peptideRT + (peptideRTWin/2))+1;

                //System.out.println("Feature:" + mass + " " + charge + " " + rtStart + " " + rtEnd);
                
                double mh = mass+PROTONMASS;
                
                String spectrumHeader = "";
                spectrumHeader += "BEGIN IONS\n";
                spectrumHeader += "TITLE=" + mh+"_"+String.format("%.2f", peptideRT) + "_" +rtStart+"_"+rtEnd+ "\n";
                spectrumHeader += "RTINSECONDS="+peptideRT+"\n";
                spectrumHeader += "PEPMASS=" +  mh + "\n";       //PEPMASS is mz, so in this case we are using MH+
                spectrumHeader += "CHARGE=1+\n";
                System.out.println("TITLE=" + mh+"_"+String.format("%.2f", peptideRT) + "_" +rtStart+"_"+rtEnd);
                HashMap<String,Object[]> uniquePeaksAt1DP = new HashMap();
                double fragRT = 0.0;
                HashMap<String,Double> uniqueApexToRT = new HashMap();
                
                String spectrum = "";
                int i = rtStart;
                while(i <= rtEnd){
                    Integer scanNum = rtIntegerToMSNScanNum.get(i);
                    i++;
                    if(scanNum!=null){
                                            
                        fragRT = Double.parseDouble(scanNumToRT1DP.get(scanNum));
                        ArrayList<Object[]> peaks = scanNumToAllDechargedData.get(scanNum);                        
                        if(peaks!= null){ 
                            //This code loops through all the peaks so we only have one entry (most abundant) per mass value at 2 decimal places
                            for(Object[] peak : peaks){
                                double ms2Mass = (double)peak[0];
                                int ms2Intensity = (int)peak[1];
                                
                                
                                //String massString = String.format("%.2f", ms2Mass); //Format to 2dp
                                
                                String massString = df.format(ms2Mass); //Format to 1dp
                                
                                
                                //Double.parseDouble(df.format(
                                
                                
                                Integer intensity = 0;
                                if(uniquePeaksAt1DP.containsKey(massString)){
                                    Integer oldIntensity = (int)uniquePeaksAt1DP.get(massString)[1];
                                    if(ms2Intensity>oldIntensity){                                              
                                        uniquePeaksAt1DP.put(massString,peak);
                                        uniqueApexToRT.put(massString, fragRT);
                                        System.out.println("\tFound higher intensity mass string: " + massString);
                                    }
                                }
                                else{
                                    uniquePeaksAt1DP.put(massString,peak);
                                    uniqueApexToRT.put(massString, fragRT);
                                    System.out.println("\tNew mass string: " + massString);
                                }  
                            }
                        }
                    }
                }
                
                for(String massString : uniquePeaksAt1DP.keySet()){
                    Object[] peak = uniquePeaksAt1DP.get(massString);
                    double ms2Mass = (double)peak[0];
                    int ms2Intensity = (int)peak[1];
                    double ms2MH = ms2Mass + PROTONMASS;    //Convert mass to MH+

                    //Here we check that the fragment mass is less than the peptide mass
                    //TODO may also want to check that charge is less than or equal, but requires a bit more coding
                    
                    //RT_WIN_FOR_MATCHING_FRAG_APEX
                    double ms2RT =uniqueApexToRT.get(massString);                     
                    double rtDelta = Math.abs(ms2RT - peptideRT);
                    
                    if(rtDelta > RT_WIN_FOR_MATCHING_FRAG_APEX){
                        System.out.println("\tReject frag RT" + ms2RT + " peptideRT: " + peptideRT);
                    }
                    
                    if(ms2MH < (mh-1) && rtDelta <= RT_WIN_FOR_MATCHING_FRAG_APEX){
                        spectrum+= ms2MH + " " +ms2Intensity+"\n";
                    }

                }
                if(!spectrum.equals("")){
                    writer.write(spectrumHeader);
                    writer.write(spectrum);
                    writer.write("END IONS\n\n");
                }
                

            }
            writer.close();
        }
        catch(IOException e){
            e.printStackTrace();
        }
        
    }
    
    
    /*
    public void XIC(String xmlFile, double rtStarts, double rtEnds, double mzStarts, double mzEnds) {
        final String xmlfile = xmlFile;
        final double frtStarts = rtStarts;
        final double frtEnds = rtEnds;
        final double fmzStarts = mzStarts;
        final double fmzEnds = mzEnds;
       
        final DefaultTableModel model = new DefaultTableModel();
        jTable2.setModel(model);
        model.addColumn("RT");
        model.addColumn("Intensity");

        final ProgressBarDialog progressBarDialog = new ProgressBarDialog(this, true);
        final Thread thread = new Thread(new Runnable(){
            @Override
            public void run(){
                //... Progress Bar ...//
                progressBarDialog.setTitle("Loading values");
                progressBarDialog.setVisible(true);
            }
        }, "ProgressBarDialog");
        thread.start();
        new Thread("LoadingThread"){
            @Override
            public void run(){                
               
                File file = new File(xmlfile);    
                if (file.exists())
                {                                
                    MzMLUnmarshaller unmarshaller = new MzMLUnmarshaller(file);

                    //... Iterate in the rt dimension ...//
                    MzMLObjectIterator<Spectrum> spectrumIterator = unmarshaller.unmarshalCollectionFromXpath("/run/spectrumList/spectrum", Spectrum.class);
                    while (spectrumIterator.hasNext()){
                        Spectrum spectrum = spectrumIterator.next();        

                        //... Identify MS1 data ...//
                        String mslevel = "";                        
                        List<CVParam> specParam = spectrum.getCvParam();      
                        for (Iterator lCVParamIterator = specParam.iterator(); lCVParamIterator.hasNext();){
                            CVParam lCVParam = (CVParam) lCVParamIterator.next();
                            if (lCVParam.getAccession().equals("MS:1000511")){
                                mslevel = lCVParam.getValue().trim();
                            }
                        }            
                        if (mslevel.equals("1")){
                            //... Get rt from spectrum ...//
                            double rt = 0.0;
                            String unitRT = "";
                            List<CVParam> scanParam = spectrum.getScanList().getScan().get(0).getCvParam();
                            for (Iterator lCVParamIterator = scanParam.iterator(); lCVParamIterator.hasNext();){
                                CVParam lCVParam = (CVParam) lCVParamIterator.next();
                                if (lCVParam.getAccession().equals("MS:1000016")){
                                    unitRT = lCVParam.getUnitAccession().trim();
                                    if (unitRT.equals("UO:0000031")) //... Validating rt unit (mins or secs) ...//
                                    {    
                                        rt = Float.parseFloat(lCVParam.getValue().trim()) * 60;
                                    }
                                    else{
                                        rt = Float.parseFloat(lCVParam.getValue().trim());
                                    }
                                }    
                            }
                            //... Get XIC across intervals ...//
                            if (rt>=frtStarts&&rt<=frtEnds){
                                Number[] mzNumbers = null;
                                Number[] intenNumbers = null;
                                //... Reading mz Values ...//
                                List<BinaryDataArray> bdal = spectrum.getBinaryDataArrayList().getBinaryDataArray();
                                for (BinaryDataArray bda:bdal){
                                    List<CVParam> cvpList = bda.getCvParam();
                                    for (CVParam cvp:cvpList){                            
                                        if(cvp.getAccession().equals("MS:1000514")){
                                            mzNumbers = bda.getBinaryDataAsNumberArray();                                
                                        }
                                        if(cvp.getAccession().equals("MS:1000515")){
                                            intenNumbers = bda.getBinaryDataAsNumberArray();                                
                                        }
                                    }
                                }                        
                                //... Generating XIC ...//
                                double intensXIC = 0.0;
                                if(mzNumbers!=null){
                                    for (int iI=0; iI<mzNumbers.length;iI++){
                                        if(mzNumbers[iI].doubleValue()>=fmzStarts&&mzNumbers[iI].doubleValue()<=fmzEnds){
                                            intensXIC+=intenNumbers[iI].doubleValue();                                    
                                        }
                                    }
                        model.insertRow(model.getRowCount(), new Object[]{
                            Float.parseFloat(String.format("%.2f", rt)),
                            Float.parseFloat(String.format("%.4f", intensXIC))});
                                }
                            }
                        }
                    }        
                    progressBarDialog.setVisible(false);
                    progressBarDialog.dispose();
                }
                else
                {
                    progressBarDialog.setVisible(false);
                    progressBarDialog.dispose();                    
                    JOptionPane.showMessageDialog(null, "File not found.", "Error", JOptionPane.ERROR_MESSAGE);
                }
            }
        }.start();
    }
*/
    
    
    
}
